Tkachenko waves in rapidly rotating Bose-Einstein condensates 
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We present a mean-field theory numerical study of Tkachenko waves of a vortex lattice in trapped 
atomic Bose-Einstein condenstates. Our results show remarkable qualitative and quantitative agree- 
ment with recent experiments at JILA. We extend our calculations beyond the conditions of the 
experiment, probing deeper into the incompressible regime where we find excellent agreement with 
analytical results. In addition, bulk excitations observed in the experiment are discussed. 
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The dynamics of quantized vorticies are important in 
a wide variety of physical phenomena from turbulence to 
neutron stars, to superconductivity and superfluidity 
The broad relevance of this problem is partly responsi- 
ble for the current interest in vortex lattices in trapped 
Bose-Einstein condensates (BECs). This system is also 
attractive because it can easily be manipulated and im- 
aged 0, IS H and because the condensates is accurately 
described within mean-field theory. Until recently much 
of the progress made in the study of vortex arrays in 
neutral superfluids was in the context of 4 He. This sys- 
tem however is notoriously difficult to manipulate and 
control. Vortex lattices in BECs have thus created an 
unprecendented opportunity to study vortex matter in 
detail. This has been recently demonstrated in spectac- 
ular experiments at JILA in which Tkachenko waves as 
well as hydrodynamic shape oscillations were directly ob- 
served y . In this paper we present a numerical study of 
these vortex matter excitations. Our approach has the 
advantage that it does not require any assumptions about 
the compressibility or homogeneity of the condensate. 

In the limiting case of an incompressible irrota- 
tional fluid, Tkachenko waves are transverse vortex- 
displacement waves traveling in the triangular lattice of 
vortex lines which constitute the ground state of the ro- 
tating superfluid In this idealized situation, the size 
of the vortex core is infinitesimal and the dynamics of the 
fluid density may be ignored. For a simple illustration, 
consider a system of n vortices each of circulation k, in an 
irrotational, incompressible fluid of density p contained 
by a vessel rotating at angular frequency f2. Let v (r) 
represent the fluid velocity at position r. If the position 
of the ith vortex is labeled by a complex number Zi, the 
free energy (/) of the array may be written as 
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Here we ignore the effects of the boundaries and the en- 
ergy associated with vortex cores pHHOI- By taking into 
account the irrotational (V x v = 0) and imcompressible 
(V • v = 0) nature of the fluid, we may derive this ex- 
pression from the simple classical relation for the energy 



which in this case is purely kinetic 
f = P 
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We can observe from Eq. that the vortex lattice, 
which represents a local minimum of this free energy, 
results from a competition between two terms: the loga- 
rithmic term which represents the intervortex repulsion, 
and the quadratic rotational term which pulls the vortices 
towards the center. The Tkachenko waves are organized 
oscillations of the vortex positions about this equilibrium. 
On the microscopic level, the vortices precess about their 
equilibrium position in elliptical orbits against the trap 
rotation with the major axis perpendicular to the wave 
propagation. On the macroscopic scale, the array is seen 
to undergo harmonic distortions which shear the lattice 
and cause the rotation of the array to alternatively slow 
down and speed up. 

Vortex arrays produced in rapidly rotating trapped 
BECs depart from this idealized case in two important 
aspects: the system is significantly compressible, and the 
underlying fluid is inhomogeneous. Due to finite com- 
pressibility, it is important to consider the dynamics of 
the underlying fluid for an accurate treatment. This 
is especially true in the light of recent interest in lat- 
tices in the mean-field quantum Hall regime 0, • To 
date, most vortex lattice studies rely on incompressible 
hydrodynamic approximations 0, 0, H IToj - The hy- 
drodynamic approximation restricts the calculation to 
longwavelength excitations and the assumption of in- 
compressibility ignores the dynamics of the fluid density. 
Baym has recently demonstrated the importance of 
compressiblity to the description of current experimental 
data. The compressibility of the lattice may be gauged 
from the ratio £ 2 /e 2 . Here £ and e are the healing length 
and intervortex distance respectively. The healing length 
£ estimates the size of the vortex core. We define these 
as £ = l/y/8nna and e 2 = 2tt/ (\/3f2) > where a is the 
scattering length of the trapped species, n is the den- 
sity at the center of the condensate and f2 is the angular 
frequency of the of the trap rotation. In this Letter we 
bypass these familiar approximations by adopting a full 
numerical approach based on mean-field theory. 

As an initial step, we must construct the undisturbed 
vortex array to a high degree of accuracy. We focus on 
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the excitations of triangular vortex arrays with hexagonal 
symmetry, such as those obtained by recent experiments. 
A pancake shaped trap geometry is employed. We note 
that owing to the centrifugal reduction in the radial trap 
frequency u> r , this is a reasonable assumption at high ro- 
tation frequencies where most of our studies take place. 
More importantly, the pancake geometry has been exper- 
imentally justified by the high quality of direct imaging 
of high [T3 and low energy lattice excitations - if the 
curvature of vortex lines were significant, such imaging 
would not be possible. According to the T — mean field 
theory, the condensate wavefunction in a frame rotating 
with angular velocity f2 satisfies the equation 
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where m is the mass of the species, uj r is the radial trap- 
ping frequency, N is the number of trapped atoms and 
L z is the angular momentum operator defined as —ih-J^. 
The coupling constant g and the external potential V(r) 
are defined as g — Amah 2 jm and V(r) = \mixi 2 r 2 respec- 
tively. The corresponding stationary equation is solved 
using a steepest descent technique starting from a trian- 
gular seed array. We define the fluctuation of the order 
parameter Sip by, 



Sip — u{r)e" 



v(ry 



(4) 



Following the Bogoliubov-de Gennes procedure, first or- 
der expansion around the ground state in terms of Sip 
then yields the eigen-equation for the exciation ampli- 
tudes w(r) and v(r) |l2(: 
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Unless otherwise stated, we express energy, length and 
time in the oscillator units hu) r , ^/h/muj r and l/uo r re- 
spectively. For the solution of the coupled eigenvalue 
equation(Eq. (J5J), we use a finite element method [l4| . 
The eigenvalue problem is characterized by three length 
scales: the Thomas- Fermi radius Rtf = V^M; the heal- 
ing length £, and the inter- vortex spacing e. Rtf specifies 
the extent of the lattice, and £ and e constrain the den- 
sity of the grid required by the calculation. Care must 
be taken to ensure the grid captures features over these 
length scales and this renders the problem very compu- 
tationally expensive. It was not unusual to obtain sparse 
matrices of the order 10 5 xl0 5 ; to solve such large eigen- 
systems, we used a variety of numerical libraries and vi- 
sualization tools Ual. 

Following Ref . [10| , we represent a Tkachcnko wave by 
a complex valued function D n rn (r, t) defined by 

D n>m (r, t) = e^[.C m (r) sin (m</> - ut - a) (6) 
+ifn. m { r ) cos (W> - ut - a)} 

The integer pair (n, m) represent the radial and angular 
order of the excitations, r = (r,tfi), a is an arbitrary phase 
and {(fn m (r), fn,m( r ))} are rea l valued functions. The 
real and imaginary part of -D„ !m , represent the displace- 
ment of a vortex at r along the x and y axis respectively. 
According to D n m , every vortex moves in an elliptical 
path. 

One of the challenges of comparing our numerical cal- 
culation to experiment or to analytic results is the proper 
identification of the excitations obtained. In our calcula- 
tion, the degree of freedom is the fluctuation of the field 




FIG. 1: Both figures depict a snapshot of the (1,0) mode in 
action, a) Simulation, b) Result from the JILA experiment. 
The curve is a fitted sine wave of wavelength 1.33Rtf 



amplitude dip, while in analytical work it is the distor- 
tion -Dn.m of the vortex array that is employed. One 
approach, involves making "movies" of each excitation 
by using Eq. (0J. In Figs. ^ and [3 we show snapshots 
from two such movies for the (1, 0) and (2, 0) modes and 
compare them to experimental data. This technique 
however, has limitations. Even with a high-resolution 
movie, the excitations are usually difficult to tell apart 
by eye, especially as £ 2 /e 2 grows. The largely trans- 
verse Tkachenko waves become increasingly corrupted by 
sound waves traveling in the underlying fluid which scat- 
ter off the vortices. A more rigorous and reliable ap- 
proach involves extracting the distortion D n ^ m from the 
fluctuation dip induced by a particular excitation. 

We therefore adopted a quantitative approach based 
upon determining the instantaneous position of each vor- 



FIG. 2: (2,0) mode in action, a) Simulation, b) Result from 
the JILA experiment. In both cases, the curves are drawn to 
guide the eye. 
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FIG. 3: Energy of the (1,0) mode vs. £ 2 /e 2 . In Ref. O the 
energy of this mode is given as 1.43ef2. 



tex. A coarse distortion function D n rn can then be con- 
structed by interpolating on the triangular lattice which 
specifies the groundstate. Along any circle centered at 
the origin, r is constant and the interpolated distor- 
tion D n ,m defines a one dimensional function d(<f>). The 
Fourier transform of e~ l ^d ((/)) then yields a peak near 




FIG. 4: Distortion of an array (£ 2 /e 2 = 0.007) in a (1,0) 
mode. The dots represent the equilibrium positions of the 
vortices, and the straight arrows represent the actual calcu- 
lated vortex displacements. The curve is a fit of a sine wave 
of wavelength 1.33-Rtf to the vortex displacements. 




FIG. 5: Distortion of an array (£ 2 /e 2 = 0.032) in a (1,0) 
mode. Compare with Fig. [I] where (£, 2 /(- 2 ) is much smaller. 
All symbols have the same meaning as in Fig. [I] 



the value m for the excitation labeled by (n,m). This 
analysis reveals one important result: the sense of pre- 
cession of each vortex in an m =/= mode may either be 
with or against the trap rotation. In these modes the 
sense of precession of a vortex is a function of its radial 
coordinate only, and along any circle all vortices precess 
in the same sense. This is a dramatic departure from 
the ideal situation of an array in a uniform, irrotational 
and incompressible fluid. In that case the vortices pre- 
cess only against the trap rotation. We note that this 
has been analytically examined |l6j . 

Our central result is illustrated by Fig. |21 where we 
make a quantitative comparision with the JILA experi- 
ment 4] and with recent analytic results 0. We focus 
on the (1,0) mode. As the parameter £ 2 /e 2 is lowered 
(abscissa), we observe a quantitative deviation from an- 
alytic results (unity line on the ordinate) which is indis- 
tinguishable from that experimentally reported |4|. In 
proportion to this deviation, we find that the eccentric- 
ity of the elliptically polarized vortex oscillations in the 
Tkachenko wave is reduced and the individual vortex mo- 
tion become almost circular. We illustrate this trend in 
Figs.E]and|S](to be compared with Fig. 3 (b) in Ref. [Io|). 
Notice that the vortex displacement vectors in Fig.[S]havc 
a more significant longitudinal projection than in Fig. 0] 
in which the ratio £ 2 /e 2 is almost an order of magnitude 
larger. In the opposite limit as £ 2 /e 2 — > 0, we obtain very 
good agreement with the analytical results. We point out 
that it is not surprising that our energies begin to exceed 
the analytic value (unity line). In the analytic treat- 
ment [Ty , the energy of this mode was calculated to first 
order in e. For small £ 2 /e 2 , we consider arrays for which 
e rs 0.2 and for which higher order terms in e could be 
important. 

In the limit of an incompressible fluid, the potential 
energy of the system is quenched and Tkachenko waves 
emerge as the only nontrivial class of excitations of a vor- 
tex array. However, they are just one of several types of 
excitations which may occur within the compressible vor- 
tex arrays realized in experiments pj . Of particular are 
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FIG. 6: a) A bulk (2,0) mode in action. The curve is 
a sine wave of wavelength 1.33Rtf predicted for the (1,0) 
Tkachenko wave, b) Phase plot of Sp Eq. 10. Notice the n 
phase change as you radially cross a node, c) Approximate 
region in the array from which atoms are removed to excite 
the (2,0) bulk or (1,0) Tkachenko modes in the JILA experi- 
ment, d) Close-up of b) on the location of a vortex. In both 
of the phase plots above, phase changes from to 2n as you 
go from dark to light regions. 

the bulk modes 0, 0, , some of which can be very 
similar in appearance to Tkachenko waves. In Ref. 
the observation of a rapid mode which is visually indistin- 
guishable from the Tkachenko (1,0) mode is reported. We 
also idenitify such a mode at exactly the experimentally 
reported value of 2.2bw p - this is just the second order 
bulk breathing or (2, 0) mode viewed in a rotating frame. 
During each half of the linear oscillatory motion imposed 



on the vortices by the breathing action, the transverse 
(Magnus) force acts in opposite directions. On a micro- 
scopic level, elliptical motion of the vortices is observed. 
Macroscopically, the array twists around the center of ro- 
tation like a torsional pendulum as the density breaths 
radially. We depict this second order 'breathing' mode 
in Fig. a) . Surprisingly, the sine wave radial distor- 
tion predicted |l(J for the (1,0) Tkachenko mode is also 
in very good agreement with this (2,0) bulk mode. To 
further illustrate this point, we make use of the complex 
density fluctuation defined for the nth excitation by [17| , 



Sp = (ip*u n - v^ip). (7) 

A casual comparision of the phase of Sp for this mode in 
Fig. b) with a schematic of the experimental probe in 
Fig. c) confirms the correct identification of the mode 
excited in the experiment. A close up of this phase plot 
around the position of a vortex is shown in Fig.[|)]d). This 
describes a clockwise elliptical oscillation of the density 
fluctuation which has a local peak at the vortex core. 
A similar motion of the vortex at this site is implicated. 
We predict the existence of third and fourth order breath- 
ing modes at 2.55w p and 3.30w p respectively. This is a 
clear example of how the presense of a transverse force 
may blur the visual distinction of the largely longitudinal 
modes from the transverse Tkachenko waves Q. We have 
observed promising candidates for the differential longi- 
tudinal waves which have been suggested previously |4j. 
In these relatively high energy excitations, vortex oscil- 
lations drive and are driven by large scale fluid density 
fluctuations. We shall report on these modes elsewhere. 
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